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As part of a project to obtain better optical response functions for nano materials and other 
systems with strong excitonic effects we here calculate the exchange-correlation (XC) potential of 
density-functional theory (DFT) at a level of approximation which corresponds to the dynamically- 
screened-exchange or GW approximation. In this process we have designed a new numerical method 
based on cubic splines which appears to be superior to other techniques previously applied to the 
"inverse engineering problem" of DFT, i.e., the problem of finding an XC potential from a known 
particle density. The potentials we obtain do not suffer from unphysical ripple and have, to within a 
reasonable accuracy, the correct asymptotic tails outside localized systems. The XC potential is an 
important ingredient in finding the particle-conserving excitation energies in atoms and molecules 
and our potentials perform better in this regard as compared to the LDA potential, potentials from 
GGA:s, and a DFT potential based on MP2 theory. 

PACS numbers: 31.15.Ew, 31.25.-v, 71.15.-ro 
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I. INTRODUCTION 

The work presented here is part of an ongoing inves- 
tigation aimed at finding more accurate and computa- 
tionally efficient ways to calculate optical spectra from 
materials with strong excitonic effects. Typical examples 
are most semiconductors, rare-gas crystals, atoms and 
molecules and nano materials. In these systems the at- 
traction of the excited electron to the hole created in the 
excitation process is very important and will strongly in- 
fluence the measured spectra both in the continuum and 
in the discrete part of the spectra. This so called particle- 
hole effect or vertex correction must be accounted for in 
the theoretical description in order to have a reasonably 
accurate, quantitative agreement with the experimental 
results. 

The traditional way to incorporate the particle-hole 
attraction within many-body perturbation theory is to 
obtain the optical spectra from solutions to the Bcthc- 
Salpeter equation for the four-point vertex function. Un- 
fortunately, this is a very demanding procedure already 
in high symmetry cases like crystalline solids and the 
method becomes very tough indeed when the symmetry 
is much lower like, e.g., at surfaces or, even worse, in 
nano systems. In solids, the two-body wave function de- 
scribing the particle and the hole is usually expanded in 
valence- and conduction-band one-electron orbitals lead- 
ing to huge secular problems. 

Nevertheless, the Bethe-Salpeter approach has pro- 
vided a large number of accurate and useful results in 
many different systems in spite of the fact that many 
approximations enter the theory like, e.g., that the ac- 
tual particle-hole interaction is most often replaced by a 
statically screened potential. During the past ten years, 
Time-Dependent Density-Functional Theory (TDDFT) 
has emerged as a competing and, perhaps, more efficient 
method for calculating optical absorption spectra includ- 
ing excitonic effects. The central quantity of this ap- 



proach is the so called exchange-correlation (XC) kernel 
the knowledge of which, unfortunately, is rather limited. 
The kernel is the functional derivative with respect to 
the density of the XC potential which is also a relatively 
unknown quantity. Due to the pioneering work of sev- 
eral groups within some of the research networks of the 
European Unioni^ it has been shown that the accurate 
results of the Bethe-Salpeter approach can be reproduced 
within the framework of TDDFT provided an appropri- 
ate XC kernel is used in the calculation. These results are 
very encouraging, the drawback being that the method 
does not provide a recipe for improving the kernel with- 
out first improving the underlying approximations within 
the Bethe-Salpeter approach. But there is a great compu- 
tational advantage of any approach based on TDDFT as 
compared to a more standard many-body approach like, 
e.g., that based on the Bethe-Salpeter equation. The for- 
mer theories are based on two-point correlation functions 
whereas the latter on four-point correlation functions. 

The simplest approximation to the XC kernel is the 
so called adiabatic local-density approximation (ALDA) 
in which the potential is taken to be the exchange- 
correlation part of the chemical potential of the electron 
gas evaluated at the instantaneous electron density. This 
simple approach is, however, known to suffer from many 
ailments. The accurate results of the Bethe-Salpeter ap- 
proach can, e.g., not be reproduced. The next level of ap- 
proximation is the exchange-only (EXO) approximation 
which, to our knowledge, has only been applied once to 
solids, 3 unfortunately, with rather poor results as far as 
the spectral properties are concerned. The approxima- 
tion seems, however, to give reasonable excitation ener- 
gies in localized systems and the total energies calculated 
from the corresponding response function are rather ac- 
curate in all systems 4^ 

Until recently, one of the drawbacks of TDDFT has 
been the lack of a systematic approach for obtain- 
ing successively better approximations to the XC ker- 
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nel. Straight-forward many-body perturbation theory 
(MBPT) can certainly be used to generate approxima- 
tions to the electronic self energy and the three-point ver- 
tex function. But the subsequent conversion into an XC 
kernel of TDDFT is not at all guaranteed to yield, e.g., 
a particle conserving density response function.—'? Parti- 
cle and current conservation are important properties of 
physical response functions which should be built-in to 
their construction. A systematic theory for improved ker- 
nels within TDDFT has recently been introduced by usi£ 
This theory automatically leads to conserving response 
functions. It is based on an adaptation to TDDFT of our 
variational approach to many-body perturbation theory^ 
A variational functional of the one-electron Green func- 
tion is constructed which yields the total energy of the 
system when evaluated at that Green function which ren- 
ders the functional stationary. Because of the stationary 
property of the functional and thus an absence of first- 
order errors, accurate energies can be obtained already at 
rather crude Green functions like, e.g., a non-interacting 
one^ The construction of the functional has two basic 
ingredients: i) a choice of basic functional expression ulti- 
mately responsible for the variational quality of the total 
functional (the size of the second-order corrections to the 
energies). So far, only two such basic expressions have 
been considered, one due to Luttinger and Ward 1 — and 
a simpler one due to Klehiji 2 - The former was shown to 
have better variational properties than the latter ji° How- 
ever, as we have shown previously— the LW functional 
leads to a response function the calculation of which is 
beyond our present day capabilities in realistic systems. 
Therefore, we will here consider only the Klein version 
of the functional, ii) The choice of $ functional which 
is also a functional of the one-electron Green function 
G. The formal significance of the ^-functional is that its 
functional derivative with respect to G yields the elec- 
tronic self-energy and its physical significance is that it 
will contain our physical intuition concerning the impor- 
tance of different physical processes. For instance, in an 
extended system all the screening diagrams should be in- 
cluded in the ^-functional and in a system with strong 
correlation effects it would be appropriate to include 
some particle-hole and particle-particle ladders. In the 
present study we will focus on the screening diagrams in 
the ^-functional and the resulting electronic self-energy 
will thus be that of the GW approximation*^ Some con- 
sideration will, however, be given also to the second-order 
exchange diagram. The variational approach to MBPT is 
converted into a density-functional theory by a restriction 
of the variational freedom for the Green function which 
is taken to be one pertinent to a non-interacting system 
moving in some local external potential V(r). Due to 
the Hohenberg-Kohn theorem applied to non-interacting 
electrons and the one-to-one correspondence between the 
applied potential and the particle density this restriction 
immediately turns the variational functional of the Green 
function into a variational functional of the density. A 
rewriting of the Klein functional in terms of the particle 



density recovers the normal form of the total energy of 
DFT— in which the so called XC energy becomes the 
$ part of the functional. Functional differentiation with 
respect to the density yields the XC potential which is 
the central object of interest in the present work. 

Worked aimed at finding the correct density-functional 
(DF) potential pertinent to different approximations 
started a long time ago. Even before the advent of DFT, 
Sharp and Horton 1 ^ proposed a method for finding a lo- 
cal potential which would accurately reproduce the total 
energies and densities of atoms within the Hartree-Fock 
approximation. This work is, in todays language, best 
referred to as the first appearance of the Exchange-Only 
(EXO) approximation or the Exact-Exchange (EXX) 
method. This method, sometimes also referred to as an 
optimized potential method (OPM) was later used by 
Talman and Shadwick^ for doing calculations on many 
atoms. In 1982, one of us made use of the Hohenberg- 
Kohn theorem to find that local potential which exactly 
reproduces the Hartree-Fock densities of several atoms.- 7 
These calculations are not equivalent to the EXO but the 
results are, numerically, very close indeed. Shortly after- 
wards, this numerical fitting procedure was generalized 
to include also all correlation effects thus producing 18 
the first exact DF XC potentials for several atoms. As 
a matter of fact, the exact XC potential for the Helium 
atom was known already at that time through the work 
of Smith, Jagannathan and Handler^ Aryasetiawan and 
Stott3>2I also found the exact XC potential of DFT for 
some smaller atoms using a presumably more accurate 
approach which had the additional advantage of offering 
insight into the so called v-representability problem of 
DFT. In solids, early progress toward an exact DF XC 
potential was made by Godby et ah 22 ' 23 who actually 
constructed their XC potential as a solution to the Lin- 
earized Sham-Schluter equation (LSS). 23 The self-energy 
of their choice was again that of the GW approximation 
(GWA). Through the formal proofs of the present work 
we now know that their potential was in fact the full RPA 
XC potential of DFT for the semiconductors they stud- 
ied. Consequently, the work of Godby et al. was almost 
identical in spirit to that of the present work, albeit in 
solids. 24 

From the middle of the eighties the number of publi- 
cations describing work aimed at finding improved XC 
potentials for DF calculations increased rapidly - both 
with regard to approximations for use in practical calcu- 
lations and with regard to a fundamental understanding 
of the behavior of the potentials in exact cases and in 
model systems. The most accurate XC potentials pro- 
duced thus far are probably those published by Umrigar 
and Gonzei 2 ^ 

In the present work we concentrate on the XC poten- 
tial at the GW level which we here prove to be identical 
to that which minimizes the standard expression for the 
total energy within the Random Phase Approximation 
(RPA). Our main motivation is an interest in density 
and current response functions beyond the EXO and the 
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XC potentials then constitutes one of the basic ingredi- 
ents. Through the work of, e.g., Petersilka et al. 26 it has, 
however, long been realized that the accuracy of the XC 
potentials is crucial for obtaining accurate excitation en- 
ergies from TDDFT. Thus, it is certainly of interest to see 
how well our GW-based potentials perform in this con- 
text. In addition, these potentials provide a good testing 
ground for our new numerical approach based on splines 
as basis functions for electronic structure calculations in 
atoms, molecules, and solids. 

The paper is organized as follows. In Sec. II we shortly 
present the formal framework based on the variational 
approach to MBPT. In Sec. Ill we present our new nu- 
merical approach and discuss its advantages and short- 
comings. Numerical results for several spherically sym- 
metric atoms are given in Sec. IV. We discuss the behav- 
ior of the GW/RPA potentials and compare their perfor- 
mance to that of potentials of other approximations, like, 
e.g., the EXX and the MP2. We also calculate particle- 
conserving excitation energies using our calculated po- 
tentials in conjunction with the approximate XC kernel 
of PGG4 Finally, in Sec. V, we draw our conclusions 
as well as advertise our forthcoming work on response 
functions. 



II. CONSERVING APPROXIMATIONS WITHIN 
TDDFT 

Physical observables of a system of interacting elec- 
trons can be calculated within MBPT, where the cen- 
tral quantity is the one-particle propagator, or the Green 
function G. The latter has a diagrammatic expansion 
in powers of the Coulomb interaction which, in extended 
systems, always must be carried to infinite order. Guided 
by physical intuition, approximations for G can be con- 
structed by including only a selected set of diagrams, 
appropriate to the system studied. The expansion of G 
can be written in terms of Dyson's equation, 

G = G H + G H SG, (1) 

where Gh is the Hartree Green function and X is the self 
energy which contains all the many-body effects above 
the Hartree level. The Hartree Green function Gh is 
the one-electron propagator of non-interacting electrons 
moving in the total potential consisting of the external 
potential w and the Hartree potential, i.e., the Coulomb 
potential from the total electronic charge density. 

Within MBPT, also the density response function \ 
has an expansion in powers of the Coulomb interaction. 
Choosing only a subset of diagrams, albeit an infinite 
subset, results in an approximate response function which 
only by pure chance will obey important conservation 
laws like, e.g., particle number, momentum and energy 
conservation. A scheme to construct approximations 
within the framework of MBPT which are conserving was 
first proposed by Kadanoff and Baynii 27 i 28 They made 



use of a functional $[G] with the property that its func- 
tional derivative with respect to G is the self energy^ 

= <5G' ® 

As a consequence, the self energy S will be a functional 
of the interacting Green function G. The response to 
external perturbations of a system treated within such 
an approximation will involve the derivative of the self 
energy with respect to G and consequently a symmet- 
ric second derivative of <E> with respect to G. It can be 
shown that this symmetry is a sufficient condition for the 
conserving properties of the resulting response function. 
Approximations generated from this scheme are called 
^-derivable. 

About the same time, Kleini 2 - constructed a variational 
functional of G composed of the ^-functional and some 
additional terms, 

iY K [G] = $[G] - Tr {GGjj 1 - 1 + ln(-G -1 )} - iU D [G}. 

Here, U Q is the classical Coulomb interaction energy be- 
tween the electrons given by U — h J nvn. When this 
Klein functional is varied with respect to G it is seen to 
be stationary when G obeys Dyson's equation, Eq. (JTJ) . 
Moreover, at the stationary point the functional takes 
the value of the ground state energy of the system. The 
functional is general and applies to any system; the refer- 
ence to the particular system is contained in the Hartree 
Green function, Gh- Other functionals of G and $ have 
been designed like, e.g., the LW functional or the ABL 
functional.— These functionals have different and gen- 
erally better variational properties as compared to the 
Klein functional and we refer the reader to Ref. for 
a more comprehensive discussion. In this work, however, 
we will focus our attention on the Klein functional. 

Starting from the Klein functional, we can restrict 
the variational freedom of the Green function to non- 
interacting ones, G s , generated by a local multiplicative 
potential, V. The Klein functional then becomes a func- 
tional of that potential. From the Hohenberg-Kohn the- 
orem there is a one-to-one mapping between the particle 
density and the potential which turns the Klein func- 
tional into a functional of the density and our theory 
into a time dependent density functional theory. The 
simplicity of a non-interacting Green function allows us 
to convert the Klein functional into the form, 

Y K [V] = -i$[G s ] + T s [n] + Jwn + U , (4) 

where T s is the kinetic energy of non-interacting electrons 
and w is the external potential. It is now clear that — i$ 
plays the role of the XC energy, E xc *£ Varying this form 
with respect to the potential V we find it to be stationary 
when the potential is given by 

V = w + V H + v xc , 
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FIG. 1: Diagrammatic representation of the "^-functional in 
the two approximations studied in this paper. 



where Vh 
obeys 



J nv is the Hartree potential and where v x 



J i X s(h 2)v xc (2) d2 = J E(2, 3)G S (3, 1)G S (1, 2) d(23). 

^ ^ J 5 ) 

This is the well known so-called linearized Sham-Schlutcr 

(LSS) equation^ 3 , which here is seen to follow from a vari- 
ational principle rather than being just the first iteration 
of an infinite number of iterations leading to the solution 
to the full Sham-Schliiter equation ! 29 ' 30 We here remark 
that there is also a self-consistency procedure involved 
in solving the LSS because the non-interacting so-called 
Kohn-Sham (KS) response function Xs as well as the self 
energy E are both expressed in orbitals obtained from 
solving a one-electron Schrodinger equation in which the 
unknown XC potential v xc is part of the local potential. 
As pointed out by Casida, it is also worth noting that 
v xc obtained in this way can be seen as the best local 
approximation to the self energy in a variational sensed 

The conserving properties, and particle conservation 
in particular, are important when calculating response 
functions from TDDFT. Within TDDFT, the interacting 
density response function \ can be shown to be given by 3 - 2 . 



X = Xs +Xs[v + fxc]X> 



(6) 



where 
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A further variation of Eq. [5] with respect to the den- 
sity gives us an equation for the two-point XC kernel, 
/ xc 4 And, because of the underlying $-derivability of the 
theory the resulting response function obeys, e.g., par- 
ticle conservation which, in the linear response regime, 
amounts to the /-sum rule. Another important property 
which also follows from the variational and conserving 
character of this theory is the well known Virial Theo- 
rem which can be used as a stringent test of the accuracy 
of the calculation of the total energy. In Appendix A the 
reader can find an explicit derivation of that theorem in 
the context of the present theory. 



A. The GW approximation 

In this work we are interested in investigating the DF 
approximation resulting from choosing $ in the GWA, 
see Fig. [Tj This approximation corresponds to a summa- 
tion of all ring diagrams, where the propagators are KS 
Green functions, 



$gw = l T r{ln(l + wG s G s )}. 
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(J) 



Inserting <j> GW into the Klein functional results in an en- 
ergy expression corresponding to the well known random 
phase approximation (RPA) for the total energy. Conse- 
quently, the total energy within the RPA is a variational 
expression and the stationary point has been proven to be 
a minimum!^ 3 - Its minimization with respect to the den- 
sity yields the ground-state energy and density as well as 
the corresponding XC potential w xc . 

Taking the functional derivative of <& in Eq. |Q we 
obtain the self energy, 



E(L2) = iG s (l,2)W(L2), 
where the effective interaction W is given by 



(8) 



(9) 



and Xs is the KS non-interacting polarization propagator. 
The LSS-equation can then be split into two terms and 
written symbolically as 



J Xs^xc = J iG s [v 



vx 



,RPA 



v\G„ 



(10) 



Keeping only the first term, the Hartree-Fock term, re- 
sults in what is known as the exact exchange (EXX) 
approximation or, sometimes, the exchange only (EXO) 
approximation and has been discussed earlier by several 
people^^ i 15 ' 16 The second term gives the correlation 
part of the potential and is expressed in terms of the 
interacting polarization propagator x RPA > which, within 
the ring approximation or the RPA, is given by 



X 



RPA 



Xs + XsVX 



RPA 



(11) 



In the following we will denote the self-consistent XC po- 
tential corresponding to the GW-level of the Klein func- 
tional by v^F A . This potential is more commonly referred 
to as the XC potential of the RPA. 



B. The second order approximation 

In the previous section a summation of ^-diagrams up 
to infinite order in the Coulomb interaction is carried out. 
A conserving approximation does, however, not require 
an infinite set of ^-diagrams. In the conserving second 
order approximation, also known as the Born approxi- 
mation, all diagrams which are at most second order in 
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the Coulomb interaction and only those diagrams are in- 
cluded in the ^-functional, see Fig [TJ Except from the 
first order Fock diagram there are then, in total, two 
more diagrams to be considered. One is the first screen- 
ing diagram, also included in the GWA. The second is the 
first vertex diagram also referred to as the second order 
exchange diagram and it is not included in the GWA. 

When considering only non-interacting KS Green func- 
tions and choosing $ in the second order approximation 
the resulting energy expression is easily seen to be equiv- 
alent to the second order energy expression obtained in 
M0llcr -Plcssct perturbation theory (MP2). This approx- 
imation is common in quantum chemistry where the in- 
serted orbitals are those of the Hartree-Fock approxima- 
tion. In our variational approach also the MP2 expres- 
sion for the total energy can be minimized with respect 
to the density thus yielding a minimizing local XC poten- 
tial. This potential was recently calculated by Jiang and 
Engcl. 34 For the purpose of testing our new numerical 
approach described in the next section, we have repeated 
their calculation. In accordance with the notation of their 
paper we will denote that correlation potential ^ P2 . 



III. NUMERICAL APPROACH 

The solution to the linearized Sham-Schliiter equation 
is complicated by the singularities of the kernel Xs giv- 
ing rise to numerical instabilities. Also, \s contains an 
infinite number of unoccupied states leading to integrals 
over the continuum. When, like here, a finite basis is 
used these integrals turn into discrete sums. Previous 
atomic and molecular calculations have lead to both un- 
physical oscillations and to an incorrect asymptotic be- 
havior of the potential, as has been discussed by many 
authorsi 34 i 35 '2£'2£'2£ In an attempt to avoid these difficul- 
ties we have designed a new basis set consisting of cubic 
splines. 



A. Cubic splines as radial basis functions 

We start by distributing a set of five nodes, not neces- 
sarily equidistant, along the radial axis, 

R = {r k : k = 0, .. . ,4;r fc < r k+1 ;r k e M}. 

From these nodes we can form a localized, piecewise 
third-order polynomial function S, in the following way: 

1. Cubic polynomials P k = a k r 3 + b k r 2 + c k r + d k , 
arc defined on the four intervals I k = 

k = 1, . . . , 4. In each of these S(r) = Pfe(r). 

2. When r < ro and r > f4, the function S 1 , is zero. 

3. The function S, is required to be continuous and to 
have a continuous first and second derivative on the 



whole real axis. This means that for k — 1, ... ,3: 

Pk{r k ) = P k+ i(r k ) 
Pk(rk) = PL + i(r k ) 
Pk(r k ) = Pk+M, 

and at the end-points: 



Pi(ro) 


= 0, 


P 4 {r 4 ) = 


P{(ro) 


-0, 


P'M) = o 


P"(r ) 


= 0, 


Pl'(r 4 ) = 



Functions designed in this manner are in the numerical 
literature called cubic splines^ From the way the spline 
above is constructed there are sixteen unknown coeffi- 
cients of the four cubic polynomials and fifteen matching 
conditions. The coefficients can then be determined up 
to a common factor. By fixing this factor the spline is 
uniquely defined on the given set R. 

For the purpose of building up a basis set, we distribute 
a set of mesh points M — {r k : k = 0, . . . , N + 3; r k < 
r k +i; r k el}, along the radial axis. From the set M we 
can define the subsets Ri — {r k : k = i, . . . , i + 4; r k € R} 
and on every subset Ri we can define a spline Si. This 
generates a basis set of TV splines with a distribution 
in space determined by the choice of mesh. This mesh 
should certainly be chosen to suit the physical problem 
at hand. One of the first computer codes for atomic 
calculations was constructed by Herman and Skillman 40 
who chose radial mesh points with a separation increas- 
ing quadratically with the distance from the nucleus - a 
so called cubic mesh. In later years an exponential mesh 
became more common but the idea is similar. It is im- 
portant to stack points close to the nucleus in order to 
account for the rapid oscillations produced by the strong 
nuclear potential while a much coarser mesh is sufficient 
in the outskirts of the atom where the wave functions 
decay exponentially. The exponential mesh is extreme 
in the sense that one obtains a very accurate descrip- 
tion close to the nucleus whereas that mesh gives a poor 
description some distance away from the atom. This is 
fine if only occupied states of single atoms are consid- 
ered. MBPT requires also a reasonable description of 
the excited one-electron orbitals for which the exponen- 
tial mesh is inadequate. It would also be less appropri- 
ate if one should like to add another atom some distance 
away to form a molecule. Thus, we have here settled on 
a cubic mesh and have chosen our mesh points according 
to r k — [h(k — 3)] 3 , k = 0, . . . , N + 3 where the 'spacing' 
h is determined by the relation r max = [h(N — 3)] 3 . This 
choice is further supported by convergence tests carried 
out by Stankovski. 41 Consequently, our entire numeri- 
cal procedure has two basic parameters, the number of 
splines N and the maximum radius r max outside which 
no physics is of interest to us. We are, e.g., not inter- 
ested in highly excited states or in scattering problems. 
We are interested in low lying excitations or in higher 
excitations only to the extent that they indirectly affect 
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the low-energy excitations. Of course, the convergence 
of our results with respect to both these numerical pa- 
rameters have been thoroughly checked. But we find it 
essential to stress that, due to the completeness of the 
splines, the results must converge toward the exact re- 
sults for the chosen physical approximation when N and 
r max are made arbitrarily large. Thus, there is no need 
to discuss the dependence of the results on the quality 
of the chosen basis set and on the choice of particular 
exponents of Gaussians or of Slater functions. 

The KS equation is a second order differential equa- 
tion the solutions of which are required to have a contin- 
uous first derivative. A potential with no discontinuities 
in the form of steps gives rise to solutions with a sec- 
ond derivative which is also continuous. Consequently, 
our cubic splines fulfill the basic requirements for radial 
basis functions and they have the additional advantage 
of constraining the potential to be continuous, which is 
a property reflecting our prejudices concerning a proper 
XC potential. 

Before ending this subsection, we would like to men- 
tion a further important numerical consequence of using 
splines as basis functions. Every single spline does only 
overlap with its three closest neighbors on both sides. 
Consequently, the matrices of the corresponding secular 
problem are band matrices for which there exist numer- 
ous efficient diagonalization algorithms. 



B. Re-expansion procedures 

The general practical procedure followed here in order 
to find the XC potential pertinent to the GW approx- 
imation is as follows. A starting potential like that of 
the simple LDA is used to generate a KS non-interacting 
Green function G s . From this Green function we easily 
obtain the KS non-interacting density response function 
Xs- The time-ordered version of Xs is 



,(r,r',w) = 



/ g (r)/*(r') 

LU 2 — (fdq — if]) 1 ' 



(12) 



where q = (k, /j,) is a particle-hole index, uj is the fre- 
quency, uj q = — efc is a particle-hole excitation energy, 
and f q is an 'excitation amplitude', i.e., a product of the 
occupied KS orbital tpk and the unoccupied KS orbital 
ip^. The fact that the KS response function Xs is diago- 
nal in this 'excitation basis' allows for a very simple and 
efficient way of solving for the full RPA density response 
function x RPA of Eq. (fTT|). 

Consequently, a substantial part of the numerical ef- 
fort must be spent on finding an efficient an accurate 
way of representing the product of two one-electron or- 
bitals. Just solving the ordinary KS equations presents 
a similar problem. At every step toward self-consistency 
the electron density, being the sum of the squares of the 
occupied orbitals, must be re-calculated. In a numerical 
approach based on localized orbitals an obvious choice of 



basis for the products of the wave functions would be the 
product of the basis functions. With N basis functions, 
this means that one would use N 2 basis functions for the 
product functions meaning, e.g., that the matrix describ- 
ing the response function Xs would be N 2 x A 2 . This is 
clearly an unnecessary effort. For instance, people using 
an approach based on LMTO:s would use products of 
LMTO:s for describing Xs but would only use a subset of 
the N 2 products of LMTO:s. The actual product basis 
functions included can, e.g., be determined numerically 
by measuring the degree of variational freedom gained 
by adding one extra product basis function. 42 Unfortu- 
nately, in this way, one may be able to reduce the num- 
ber basis functions for the products by a factor of two 
or three which is not a very large gain if N is large. In 
the case of methods based on plane waves the situation 
is slightly better. If the accuracy of the wave functions 
are considered enough by including all plane waves up to 
a chosen momentum cut-off, quantities like the density 
and the excitation functions will contain, on the average, 
eight times as many plane waves as the wave functions. 

In our case, using the N cubic splines as basis func- 
tions, our excitation functions would be a sum of poly- 
nomials of degree six. But if it is sufficient to describe the 
wave functions in terms of N polynomials of the third de- 
gree, intuition suggests that the same should be true also 
for the density and the excitation functions. We thus re- 
expand products of wave functions in terms of the same 
set of cubic splines as used for the wave functions. The 
accuracy of this intuition has to be verified numerically. 
We have found that, in the case of a fixed atomic po- 
tential - i.e. no charge density or self-consistency is in- 
volved - the accuracy of the KS eigenvalues increases by 
an order of magnitude and the relative error decreases to 
10~ 5 when increasing the number of splines from twenty 
to thirty. However, including also errors arising from the 
re-expansion of the self-consistent density in terms of the 
same number of cubic splines as used for the wave func- 
tions this error increases by a factor of three at thirty 
splines. By a rather minor increase in the number of 
splines - as compared to eight times the number of plane 
waves or the number of localized basis functions squared, 
we regain the full accuracy (10 -5 ) at just below forty cu- 
bic splines. From this we conclude that the cubic splines 
constitute a superior basis set for many-body calculations 
involving two-point correlation functions. 

The so obtained functions ^ RPA and G s are finally 
used to calculate the screened interaction W and the 
self-energy £ of the GWA. And then the LSS, Eq. © 
is solved for the XC potential v xc by expanding it too in 
our cubic splines and inverting Xs expressed as a matrix 
in cubic splines. As discussed above this matrix is sin- 
gular. The physical reason for this is that the response 
of a constant potential is zero. Inverting Xs is thus not 
a unique operation. This difficulty can be circumvented 
by adding the constraint: 



lim 



: (r) = 0. 



(13) 
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Mathematically this means to invert the matrix Xs in the 
subspace orthogonal to that defined by the zero eigen- 
value of Xs- 

In all our calculations the EXX potential was first cal- 
culated and then used as a starting guess. The conver- 
gence criterion in our calculations was set to \n^ k '(r) — 
n (fc-i)( r )| < io~ 5 , and was reached in a few iterations 
for all spherical atoms up to Ar. 



IV. RESULTS FOR SPHERICAL ATOMS 

In the present paper we present results for the spherical 
atoms He, Be, Ne, Mg and Ar, see Figs. [2]|4] In the cases 
of He, Be and Ne we compare our results to existing exact 
density functional potentials^ 



0.05 
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A. The RPA correlation potential 

The correlation potentials in this work are defined as 
the difference between the self-consistent XC potentials 
calculated within the RPA, or MP2, and that of the EXX, 



v c (r) = u xc [n xc ](r) - u x [n x ](r) 



(14) 



Note that the correlation potential is defined as the differ- 
ence between potentials calculated at two different densi- 
ties. For a more detailed discussion of this point we refer 
to Sec. llVCl 

In Fig. [5] the RPA and MP2 correlation potentials for 
He, Be and Ne are presented and compared to the exact 
correlation potentials^ The characteristic shell oscilla- 
tions inherent in the exact potential can not be repro- 
duced by potentials depending explicitly on the density. 
Indeed, both LDA- and gradient corrected potentials lack 
the correct shell oscillations (see e.g. Ref. [H). The po- 
tential d^ pa , however, which is an implicit density func- 
tional through the dependence on the KS orbitals and 
eigenvalues, can be seen to reproduce these oscillations 
very well and this is also the case for the potential of 
MP2. Except at the origin, the amplitudes of the os- 
cillations of v^ PA are much closer to those of the exact 
potential as compared to the case of v^ P2 . At the origin 
the RPA potential appears to deviate the most from the 
exact potential by being too attractive for all atoms. It 
should then be remembered that the values of the po- 
tential in this region is expected to be of less importance 
due to the strong singular Coulomb potential from the 
nucleus. It is tempting to interpret the superiority of the 
RPA correlation potential to that of the MP2 in the outer 
parts of the atoms as an increased importance of screen- 
ing in this region. Similarly one might guess that short 
range correlations are more important in the interior of 
the atoms leading to a better performance of the MP2 
potential in this region. 

In order to demonstrate the effect of correlations on the 
total XC potential we compare the RPA version of v xc 
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FIG. 2: Self-consistent correlation potentials for He, Be, and 
Ne. The RPA-potential is compared to the MP2-potential 
and to the exact potential. For Be no self-consistent MP2- 
potential can be obtained. Instead, the potential was evalu- 
ated at the EXX density. For a further discussion see Ref. 

M 
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to the potential of u x of EXX in Fig. [3] (for Ne and Ar). 
As expected, the shell oscillations are damped by cor- 
relation. Due to the good performance of our numerical 
method the 1/r asymptotic behavior can be observed to a 
large radius without any fitting procedure. According to 
Niquct et alM the RPA correlation potential should have 
a polarization tail of the form — a PPA /2r 4 , where ajq is 
the static polarizability of the corresponding atom. This 
effect can, however, only be observed at a considerable 
distance from the atom, « 20 bohr radii. An accurate de- 
scription of the potential at such large distances requires 
many more basis functions as compared to our standard 
calculations. Increasing the accuracy in the asymptotic 
region, a tendency toward such a tail was, however, ob- 
served. 

Finally, in Fig. [| the RPA and MP 2 correlation po- 
tentials of Mg and Ar are displayed. The presence of a 
third s-shell, as well as a second p-shell for Ar, can be 
observed. 



B. Ionization potentials 

Within exact DFT the highest occupied orbital eigen- 
value equals the negative of the ionization potential^ 
A further test of the quality of the RPA XC potential is 
thus a comparison between the eigenenergy of the HOMO 
and the experimental ionization potential. Results are 
presented in Table |U In this table we also present the 
corresponding results obtained from the EXX and from 
MP2. The latter approximations also give the proper 
— 1/r tail of the potential which is very important for 
obtaining a reasonable ionization potential. Therefore, 
we have not considered it worthwhile to include results 
from other local potentials like, e.g., those of the LDA 
or different GGA:s. As can be seen, the ionization po- 
tentials of the RPA potential are in very good agreement 



with experiment and represent a significant improvement 
on the EXX and also on MP2. Consequently, in this re- 
gard, the i^, PA is the best performing potential presently 
known. 



C. Total energies and role of self consistency 

Using the self-consistent KS orbitals and eigenvalues 
the total energy was calculated from Eq. The results 
for different atoms and are presented in Table HT1 We see 
that the correlation energy is overestimated by almost 
a factor of two for the small atoms, leading to a too 
low total energy. This tendency within the RPA has 
previously been pointed out by several workers! 10 ' 45 The 
MP2 energy functional is seen to perform much better in 
this regard. 
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TABLE I: Ionization potentials for some atoms. The RPA 
potential produces very accurate ionization potentials com- 
pared to other potentials with the correct 1/r decay. Values 
are in Hartrees. 



Atom 


RPA 


MP2 


EXX 


Exp 6 


He 


0.902 


0.893 


0.918 


0.903 


Be 


0.354 


0.357 a 


0.309 


0.343 


Ne 


0.796 


0.657 


0.851 


0.792 


Mg 


0.297 


0.303 


0.253 


0.281 


Ar 


0.590 


0.560 


0.591 


0.579 


MAE C 


0.009 


0.040 


0.030 





"Calculated using the EXX density. 
'Experimental data taken from Ref. |34| 
c Mean average error. 

TABLE II: Total ground-state energies calculated from the 
self-consistent density. The RPA energy functional gives too 
large correlation energies, whereas the MP2 functional per- 
forms much better. 



Atom 


RPA 


MP2° 


EXX 


Exp" 


He 


2.945 


2.909 


2.862 


2.904 


Bo 


14.754 


14.697 


14.572 


14.667 


Ne 


129.143 


129.027 


128.545 


128.937 


Mg 


200.296 


200.129 


199.612 


200.059 


Ar 


527.908 


527.661 


526.650 


527.604 



"Experimental data and MP2 values are taken from Ref. I 



We have found the Virial Theorem (VT) to be a conve- 
nient test of the accuracy of our calculations. Therefore, 
in Appendix A, this theorem has been proven by us in 
the case of any conserving density functional approxima- 
tion. As a matter of fact, in our calculations, the VT is 
obeyed to within six significant digits. This is a very sat- 
isfactory test on the overall accuracy of our calculations. 
It is important to point out that the VT only holds if the 
self-consistent orbitals and one-electron eigenvalues are 
used in the evaluation of the energies. 

In order to investigate the variational properties of 
the total energy of the RPA as a function of the den- 
sity we have evaluated this functional using orbitals and 
one-electron eigenvalues also from the LDA, MP2 and 
the EXX. The results are illustrated in Fig. O In all 
cases we find a higher energy compared to the fully self- 
consistent RPA result, in accordance with our previously 
proven theorem that the RPA density functional has a 
minimum. 33 We also studied the variation of the poten- 
tial with respect to the density used for its evaluation. 
Only small changes in the potential was observed. Com- 
paring to the results of Jiang and Engel 3 ^ we conclude 
that the RPA potential is more stable with respect to 
variations in the density. 

D. Two-electron excitation energies 

The particle-conserving or two-electron excitation en- 
ergies can be obtained from the full density response 



function \ of the system. In finite systems, there are 
at least a few such excitation energies below the contin- 
uum edge which then show up as poles of %. Within 
exact TDDFT these poles are zeros of the expression 
Xj 1 — v — / xc , where v is the bare Coulomb interaction, 
/ xc is XC kernel, and \s is the non-interacting KS re- 
sponse function with poles at the differences between the 
DF eigenvalues. In a one-pole approximation, the two- 
particle excitation energies are seen to be the difference 
between an occupied and an unoccupied DF eigenvalue 
corrected by some matrix elements with respect to DF or- 
bitals of the Coulomb interaction and the XC kernel, / xc . 
It has been shown previously by Petersilka et al££- that 
the latter matrix elements have a much smaller influence 
or effect on the calculated excitation energies than the 
eigenvalue difference obtained by using different approx- 
imations to the XC potential. Of course, the presence 
of a sum over all poles will affect the actual zeros of the 
denominator of x, i.e., the excitation energies. If, how- 
ever, these zeros are well separated, as is often the case in 
the discrete part of the spectrum, also this effect is much 
smaller than the eigenvalue differences produced by dif- 
ferent XC potentials. As a consequence, an accurate XC 
potential is of vital importance for obtaining accurate 
two-particle excitation energies. 

In the Tables IIIHlVl the KS eigenvalue differences in 
RPA are compared to the exact eigenvalue differences ob- 
tained from Ref. and to the MP2 and EXX eigenvalue 
differences. The results show that the mean average error 
is significantly reduced as compared to the EXX approx- 
imation. For He and Be the MP2 approximation also 
improves the EXX values. 

The magnitude of the Hartree contribution and the / xc 
contribution to the true excitation energies are presented 
in Tables [VI] and IVIII We see that the KS eigenvalue 
differences are already very close to the true excitation 




FIG. 5: The total energy for He and Ne in the RPA cal- 
culated with orbitals and one-electron eigenvalues from the 
LDA, EXX, MP2 and the RPA. The RPA orbitals are seen 
to give the lowest energy confirming that the RPA functional 
has a minimum. 
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TABLE III: KS eigenvalue differences for He. The values of 
the RPA and MP2 improves significantly on the EXX values. 



Transition 


RPA 


Ivlrz 


£jAA 


H;xact 


ls^2s 


0.744 


0.736 


0.760 


0.746 


Is— >3s 


0.844 


0.836 


0.861 


0.839 


ls^2p 


0.775 


0.768 


0.791 


0.777 


Is— >3p 


0.855 


0.846 


0.871 


0.848 


MAE 


0.004 


0.006 


0.018 




"Taken from Ref. [H 








TABLE IV: KS eigenvalue differences for Be. The 


errors in 


the EXX ei; 


;envalue differences are reduces by a factor of two 


in both the RPA and MP2. 


Note that the MP2 values are 


calculated from the EXX density. 






Transition 


RPA 


MP2 


EXX 


Exact" 


2s^3s 


0.254 


0.253 


0.214 


0.244 


2s^2p 


0.131 


0.128 


0.131 


0.133 


2s^3p 


0.276 


0.276 


0.234 


0.269 


2s^4p 


0.332 


0.330 


0.297 


0.305 


2s^3d 


0.292 


0.292 


0.241 


0.283 


MAE 


0.011 


0.011 


0.023 




"Taken from Ref. [H 








TABLE V: KS eigenvalue differences for Ne. The RPA results 


are almost 


an order of mac 


;nitude 


better than those of the 


EXX which 


are better than those of MP2. 




Transition 


RPA 


MP2 


EXX 


Exact" 


Is— >3s 


30.639 


30.591 


30.628 


30.633 


ls^3p 


30.715 


30.652 


30.706 


30.706 


ls^3d 


30.773 


30.699 


30.766 


30.759 


2s-+3s 


1.462 


1.336 


1.526 


1.469 


2s^3p 


1.538 


1.398 


1.604 


1.542 


2s^3d 


1.595 


1.444 


1.664 


1.595 


2p^3s 


0.607 


0.492 


0.659 


0.612 


2p— >3p 


0.683 


0.553 


0.737 


0.684 


2p^3d 


0.740 


0.600 


0.797 


0.738 


MAE 


0.007 


0.111 


0.046 




"Taken from Ref. [H 



energies. Including the Hartree contribution, i.e., eval- 
uating the poles of the RPA response function, drives 
the KS values further away from the true excitation en- 
ergies. Including the / xc part cancels this error and in 
most cases we come back to a value of the same quality 
as the KS eigenvalue differences. For some transitions, 
however, e.g., the 2s— >2p transition in Be, it is necessary 
to include both the Hartree and the / xc contributions. 



V. CONCLUSIONS AND DISCUSSION 

In the present work we have calculated that local po- 
tential which through a one-particle Schrodinger equa- 
tion generates the orbitals which minimize the expression 
for the total energy within the RPA. The systems studied 
are spherical atoms. We have found that the correlation 



part of this potential - defined as the total RPA potential 
minus that corresponding to the EXX - has the effect of 
softening the shell structure produced by the potential of 
the EXX. This shell structure consist of rapid oscillations 
of the XC potential between the atomic shells. 

For several spherical atoms, we have compared the spa- 
tial dependence of our RPA potential to that of the exact 
DF potential defined to be the potential which, through 
the Kohn-Sham procedure, yields the exact electron den- 
sity of the many-body system. Such exact potentials ex- 
ist in the cases of the He, Be, and Ne atoms. We have 
found that the RPA potentials are closer to the exact 
DF potentials than the corresponding local potentials of 
MP2 theory and much closer than the potentials of the 
LDA or of any of the GGA:s. 

It is relatively difficult to judge the quality of a given 
local potential from a study of its spatial dependence. 
Within DFT, the highest occupied eigenvalue ought to 
equal the negative of the ionization potential. And we 
have found that our RPA potentials perform very well 
in this regard, much better than any traditional poten- 
tial but also better than the local potential from MP2 
theory. Our average error in the obtained ionization po- 
tentials are ~0.24 eV compared to —1.1 eV in the case of 
MP 2 theory and -0.8 eV within the EXX. We interpret 
this result as being due to the GWA providing a better 
description of correlation effects among the more loosely 
bound valence electrons being ionized as compared to the 
case of MP2. 

Another measure of the quality of the calculated RPA 
potentials can be obtained from a study of two-electron 
excitation energies. Within TDDFT, and using a single- 
pole approximation, these can be obtained as differences 
between DF eigenvalues corrected by relatively small ma- 
trix elements of the bare Coulomb interaction and the XC 
kernel. The latter corrections are relatively insensitive to 
the orbitals used in their evaluation and on the choice 
of XC kernel. Due to the accurate DF potential of the 
GWA, i.e., the RPA potential, the corresponding eigen- 
values are close to exact DF eigenvalues and the resulting 
two-electron excitation energies are also very accurate. 

In the present work, we have shown that the tradi- 
tional RPA follows from a special choice of variational 
expression for the total energy involving the GWA for 
the functional $ from which the electronic self-energy is 
obtained as a functional derivative with respect to the 
Green function. This guarantees that the density re- 
sponse function obtained by perturbing the system by an 
external potential will be conserving meaning, e.g., that 
it will obey the f-sum rule. Another consequence is the 
fact that the ground-state energies obey the Virial Theo- 
rem which is here proven explicitly for any approximation 
within TDDFT obtained from the variational approach 
to MBPT. 

The demonstrated high quality of the XC potential 
within the RPA, particularly with regard to energy dif- 
ferences, induces strong hopes that the density response 
function obtained by perturbing the system will be very 
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TABLE VI: Excitation energies for He. In the first three 
columns the difference between the poles of the response func- 
tion, in three different approximations, and the experimental 
values are presented. Note that the poles of Xs are just the 
KS eigenvalue differences. The last column gives the experi- 
mental values taken from Ref. |26| . Already the KS eigenvalue 
differences are close to the experimental results. 



Transition 


(x.r 1 


(x r PA) -i 




Exp. 


ls^2s 


-0.014 


+0.010 


+0.006 


0.758 


ls-^3s 


+0.002 


+0.013 


+0.008 


0.843 


ls^2p 


-0.005 


+0.005 


+0.002 


0.780 


Is— >3p 


+0.006 


+0.012 


+0.008 


0.849 


MAE 


0.007 


0.010 


0.006 




TABLE VII: 


Excitation energies for Be. 


The columns presents 


the same quantities as in Table IVT1 The 2s^2p KS transition 


is seen to be 


responsible for the large MAE in the first column 


and needs to be corrected by the Hart: 


ree and the / xc 


term. 


Transition 


(X,)- 1 


(x rp A) -i 


(x pg G) -i 


Exp." 


2s^3s 


+0.005 


+0.017 


+0.010 


0.249 


2s^4s 


+0.013 


+0.027 


+0.014 


0.297 


2s^2p 


-0.063 


+0.024 


-0.002 


0.194 


2s^3p 


+0.002 


+0.012 


+0.007 


0.274 


MAE 


0.020 


0.020 


0.008 





"Taken from Ref. 26 



accurate indeed. But this will be the topic of a future 
publication. 

In the present paper we have introduced a novel way 
of doing electronic structure calculations based on cubic 
splines as radial basis functions. The original motivation 
for the introduction of this somewhat unusual basis set 
was the desire to circumvent numerical difficulties associ- 
ated with the known singularities of the Kohn-Sham non- 
interacting density response function. The latter gives no 
response to a constant potential (long-wave-length limit) 
and a very small response to a very rapidly varying po- 
tential (limit of short wave length). To judge from the 
high accuracy of our results the splines appear to be ide- 
ally suited to deal with the latter problem. 

We also want to stress two more advantages of our 
numerical technique based on the cubic splines, i) the 
proven accuracy of the re-expansion of a product of two 
splines in terms of splines guarantees that the matrices 
corresponding to two-particle propagators like, e.g., re- 
sponse functions, are of the same sizes as one-particle 
propagators like, e.g., Green functions or - for that matter 
- wave functions. This property is a clear advantage over 
more standard basis sets consisting of, e.g, plane waves 
or LMTO:s. ii) The secular problem based on splines 
requires the handling of sparse matrices for which there 
exist efficient standard computer codes. Our nice expe- 
rience from using the splines on atomic problems sug- 
gests that we ought to implement similar methods also 
in molecules and solids. 

Work to apply the ideas introduced in the present pa- 



per also to calculate density response functions and phys- 
ical properties like, e.g., polarizabilities, are in progress. 
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APPENDIX A: THE VIRIAL THEOREM 

We will prove that the Virial Theorem holds for any 
conserving approximation within DFT generated from 
the Klein functional. In the full many body case the 
Virial Theorem has already been proved in a conserving 
approximation^ 

Consider the Klein functional as a functional of the 
density n: 



Y K [n] = -i$[G s ]+T s [n} + J 



wn + U . 



(Al) 



Keeping the normalization, the density is scaled with re- 
spect to the spatial coordinates, 

n x = X 3 n(Xr). (A2) 

Due to the stationary property of Yk we have that 

'dY K \n x ] 



dX 



0. 



(A3) 



i=i 



Let us see how each term in the Klein functional scales 
when we scale the density as in Eq. (|A2|) . The Hartree 
term scales linearly (y = 1 jr) 



U x = 1/2 J d 3 rd 3 r'v(r ~r')X 3 n(Xr)X 3 n(Xr') 

= 1/2 J d 3 rd 3 r'v{r - r')n{r)n{r') 
= XU , 



(A4) 



so 



dUl 
dX 



= u n 



1=1 



The kinetic energy of independent particles is a sum of 
the occupied one-particle kinetic energies, and is thus an 
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implicit functional of the density. With the density scaled 
as in Eq. (|A2[) the orbitals scale like 



^ A = A 3 /V(Ar). 



(A5) 



Inserting the scaled orbitals in the expression for the ki- 
netic energy we find 



1 OCC p 

-i^/ d\X 3 / 2 ^{Xr)V 2 X 3 / 2 ^{Xr) 

i 

OCC p 

--J2J d\X 2 ^(r)V 2 ^{r) 



AT. 



(A6) 



and 



dX 



2T, 



1=1 



The external potential energy W scales as 



d 3 rX 3 n(Xr)w(r) 

d 3 rn(r)w(r / A). 

If the external potential is Coulombic we have 
(dW v 



(A7) 



\ dX 



W. 



1=1 



The XC energy E xc is not an explicit functional of n but 
of G s . To see how E xc scales we must first determine 
how G s scales, that is, to find that Green function G l s , 
which corresponds to n l for every /. Taking into account 
that the orbitals scale as in Eq. (| A5[) we have 

-ivV(Ar) = -^V 2 A ^(Ar). 

From the KS equation, 

{-^V 2 + V(r)}^ fe (r) = s k <p k (r), 



we see that 
A 2 



A 2 



2 V 2 ^(Ar) = --{e k - V{Xv)}^ k {Xv). 

Thus, the equation that yields the scaled orbitals is 

{ _I V 2 + AV(Ar)K(Ar) = A 2 e^(Ar). 

Consequently, the Green function scales as 

G*(r,r',oj) = XG{Xr,Xr',uj/X 2 ). 
A general <I>-diagram of order n can be written 
1 



(A8) 



(A9) 



2n 



TrE n [G S ]G S , 



where £„ is a skeleton diagram of order n. In order n 
there are thus n interaction lines, 2n coordinates and 2n 
propagators. Inserting the scaled Green function we find 
that <&„ scales as 



r 2 



-Tr T, n [G l s ]G l s . 



2n 



(A10) 



A factor of l 2n comes from the number of propagators. 
The variable substitution gives us a factor of l~ 6n and 
from the interaction lines an additional factor of l n is 
obtained. There are n + 1 w-integrations giving a factor 
of Z 2 (™ +1 ) . In total we obtain a factor of l 2 ~ n . Taking the 
Z-derivative we obtain 



Summing over n yields 



(2-n)$. 



d¥ 
~dX 



1=1 



n v 

£(2 -n)#„ 

n 

Th 

2tE xc [G s ]-iU xc [G s 



(All) 



where 



U xc {G s ] = - l -J2 r bZ n [G s ]G s . 



Summing all the terms we have proved the Virial Theo- 
rem 



where 



= 2T S + U n + W + 2E XC - U xc 
= 2T S + 2T XC + Un + W + U XC , 



U x 



(A12) 
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